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ABSTRACT 


This  paper  describes  and  derives  a  numerical  integra¬ 
tion  technique  presently  being  used  at  NRL  for  problems  in 
structural  dynamics.  The  method  is  capable  of  computing 
damped  and  undamped  shock  spectra,  Fourier  sine,  cosine, 
and  regular  Fourier  integrals  (frequency  response  from 
transient  response),  and  the  inversion  of  the  Fourier  sine 
and  cosine  integrals  (transient  response  from  frequency 
response).  It  can  also  be  used  to  calculate  Fourier  series 
coefficients,  and  for  the  numerical  solution  of  nonlinear 
equations.  Several  examples  are  worked  out  in  detail  and 
some  others  calculated  byNAREC  (NRL’s  digital  computer) 
are  shown  to  present  some  idea  of  the  precision  of  the 
method.  This  report  does  not  assume  a  sophisticated 
mathematical  background  and  uses  only  those  techniques 
which  are  available  to  undergraduate  students  in  engineering. 
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A  NUMERICAL  PROCEDURE  FOR  SHOCK  AND  FOURIER  ANALYSIS 


INTRODUCTION 

In  the  calculation  of  the  response  of  linear  elastic  structures  to  transient  forces  and 
foundation  motions  it  is  often  necessary  or  convenient  to  compute  the  response  of  a  set 
of  single  degree  of  freedom  oscillators  to  the  impressed  transient.  Normal  mode  theory, 
for  example,  defines  (1-3)  a  set  of  single  degree  of  freedom  systems  whose  superposed 
responses  become  the  structural  response.  It  will  be  shown  in  a  later  section  that  the 
solution  of  certain  Duhainel  integrals  (dealing  with  the  undamped  single  degree  of  freedom 
mechanical  system)  leads  to  a  means  of  evaluating  Fourier  series  coefficients  and  to  the 
representation  of  the  Fourier  transform. 

Quite  often  only  a  graphical  record  of  the  impressed  transient  is  known  or  a  Fourier 
spectrum  of  the  solution  of  the  differential  equations  is  known  and  the  problem  is  to 
calculate  the  structural  response  (4,5).  This  is  usually  accomplished  by  electrical  analog 
techniques,  some  general  numerical  integration  method  (6,7),  or  by  graphical  or  semi- 
graphical  solutions  (1,8,9).  It  is  quite  legitimate  therefore  to  ask  “Why  present  another 
method?”  The  answer  would  be  that  although  many  procedures  are  available,  most  are 
designed  to  show  how  a  special  technique,  an  iteration,  an  integration  by  replacement  of 
the  original  differential  equation  by  one  of  finite  differences,  a  Maclaurin  series  expansion 
with  associated  assumptions  in  regard  to  the  behavior  of  the  derivative  of  the  acceleration, 
or  an  analogy  is  used  to  solve  the  problem.  It  was  thought  that  the  individuals  involved  in 
solving  these  problems  might  profit  both  timewise  in  the  calculation  of  the  response,  and 
by  an  increased  understanding  of  the  problem,  if  a  simple,  powerful,  yet  precise  technique 
which  closely  followed  an  easily  understood  physical  argument  was  presented. 

In  this  method  the  differential  equation  of  motion  for  the  oscillator  is  not  replaced  in 
the  classical  sense  by  one  of  finite  differences,  nor  is  any  iteration  required  in  the  linear 
case  to  obtain  a  solution,  even  though  damping  is  present.  The  technique  is  simply  a 
logical  extension  of  ones  previously  presented  (1)  and  has  been  used  extensively  at  the 
U.S.  Naval  Research  Laboratory  since  the  author  proposed  it  in  1956. 

This  report  deals  primarily  with  the  undamped  response  of  the  oscillator  to  a  founda¬ 
tion  motion  when  the  velocity  of  the  base  is  known.  The  numerical  integration  method  is 
also  presented  with  an  example  for  this  case  with  linear  damping,  and  the  equations  for 
the  damped  or  undamped  response  to  applied  forces  or  foundation  accelerations  are  also 
presented.  The  nonlinear  cases  are  only  briefly  mentioned,  as  they  will  be  the  subject  of 
future  reports.  Some  Fourier  series  coefficients  are  calculated  and  two  inverse  Fourier 
transform  problems  are  worked  out. 


SYMBOLS 

A  dot  over  a  symbol  indicates  differentiation  with  respect  to  time. 
An  nth  Fourier  cosine  coefficient 

An  nth  Fourier  cosine  coefficient  for  half  range  expansion 
B„  nth  Fourier  sine  coefficient 

Bn  nth  Fourier  sine  coefficient  for  half  range  expansion 
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nonlinear  dissipative  damping  term 
D(x)  with  linear  term  removed 
force  as  a  time  function 
force  at  time  denoted  by  subscript 
Fourier  cosine  transform  value  at  oj. 

Fourier  sine  transform  value  at  a. 
linear  spring  constant 
a  dummy  variable 
mass 

total  number  of  data  points 
nonlinear  resistance  function 

nonlinear  resistance  function  with  linear  term  absent 
first  forward  difference  at  n 
second  forward  difference  at  n 

Sn-S„-i 

time 

time  at  end  of  input 
step  change  of  velocity  of  value  v 
relative  displacement 
absolute  displacement 
foundation  displacement 
an  angle 

linear  damping  term  in  cfc 

deviation  or  error  at  i 

frequency,  cycles/time 

a  function  of  time 

right-hand  limit  of  f(t) 

left-hand  limit  of  f(t) 

functions  of  r,  s,  etc. 

an  increment  of  time 

damped  natural  frequency,  rad/time 

time 

ratio  of  damping  to  critical  linear  damping 
negative  of  average  value  of  a  function 
3.14159  - 


0  an  angle 

K  a  frequency,  rad/time 

oj  undamped  natural  frequency  of  linear  system,  rad/time 
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THE  SIMPLE  OSCILLATOR 

The  simple  oscillator  of  Fig.  1  has  the  following  general  differential  equation  of 
motion: 


MY  +  D(X)  +  R(X)  a  F(t)  (1) 

where  X  -  (Y  -  Z) ,  D(X)  is  a  dissipative  term,  and  R(X)  is  a  restoring  force  term. 
Equation  (1)  may  be  written  in  the  form 

MY  t  C(Y  -  Z)  +  K(Y  -  Z)  +  D(X)  +  R(X)  =  F(t)  (2) 

where  6(X)  and  R(X)  are  the  nonlinear  terms  of  the  dissipative  and  restoring  forces.  If 
the  system  is  linear,  Eq.  (2)  reduces  to 

MY  +  C(Y  -  Z)  +  K(Y  -  Z)  =  F(t).  (3) 

Assuming  that  c  <  2 /km  (i.e.,  C  <  critical),  and  using  the  notation 

w2  =  K/M,  a  =  ,  P  =  01  /l  -  a2 

2Ma) 

Eq.  (3)  may  be  written  as 

Y  +  2acuY  +  a.2  Y  =  2  aw  Z  +  a>2  Z  +  .  (4) 

There  are  two  cases  that  will  be  the  subject  of  this  report:  there  is  an  applied  force,  but 
no  foundation  motion;  and  there  is  foundation  motion,  but  no  applied  force.  Combinations 
of  solution  for  the  linear  system  can  of  course  be  determined  by  superposition. 


“ I” 

0(X| 

Fig.  1  -  The  simple  oscillator 
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For  relative  motion  (X)  in  the  case  of  foundation  motion  with  no  applied  force,  Eq.  (4) 
reduces  to 

X  +  2awX  +  w2  X  =  -z  (5a) 

for  the  damped  case,  and 


X  +  u1  x  =  -z 


(5b) 
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(or  the  undamped  case.  For  applied  force  and  no  foundation  motion  they  become 

2 

/V  T  O  ^  /%  T 

and 


X  +  ^ 


X  +  J  X  = 


F(t) 

M 


(5  c) 
(5d) 


since  X  =  Y  -  Z,  and  Z  ■  0. 


DERIVATION  OF  METHOD 


Undamped  Response 

In  order  to  shorten  this  paper  only  the  numerical  integration  equations  for  a  linear 
oscillator  responding  to  a  foundation  motion  will  be  derived  in  detail.  The  problem 
considered  is  to  obtain  the  relative  response  if  a  record  of  the  velocity  of  the  foundation 
is  known  as  a  function  of  time. 


The  solution  of  Eq.  (5b)  and  its  derivative  can  be  shown  to  be  (1) 


t 


*0  1 

X  =  cos  cot  +  —  sin  wt  -  — 

°  CO  (O 

J  Z(T)  S  in  CD  (t 

-  T)  dT 

(6a) 

and 

,  1 

*0 

X  v  .  x°  „  I  f 

—  =  -X.  sin  cot  +  ■ —  cos  a»t  -  —  I 

CD  O  CO  CO  1 

2(T)  cos  co  (t 

-  T)  dT. 

(6b) 

'o 


These  are  linear  equations,  so  if  the  solution  were  known  at  t  =  t,  ,  the  values  of  displace¬ 
ment  and  relative  velocity  could  be  considered  to  be  “initial  values”  for  time  redefined  to 
begin  there.  Therefore  if  the  Duhamel  integrals  could  be  integrated  with  good  precision 
over  a  time  interval,  Att  =  h,  say,  then  these  new  values  of  x  and  x/wat  t  =  tt  +  Atj  could 
be  used  as  “initial  values”  to  continue  the  problem. 


Consider  a  portion  of  a  velocity  record  as  shown  in  Fig.  2.  For  convenience  the 
record  is  divided  in  a  number  of  regions  by  equally  spaced  time  Increments,  h  =  At.* 
Assume  that  X  and  x/ware  known  for  t  =  nh,  i.e.,  Xn  and  Xn/w.  The  problem  resolves 
itself  to  finding  the  values  of  t  =  (n  +  l)h.  If  xn  and  xn/o>are  considered  to  be  “initial 
values*  and  time  redefined  as  starting  at  this  point  only,  the  integrals  need  be  evaluated 
at  the  new  t  =  h.  Suppose  it  was  desired  to  use  the  average  foundation  acceleration  over 
the  Increment  as  the  approximation  to  the  true  foundation  acceleration;  then 


Average  foundation 
accel era  t ion 


Z(t)  dt  = 


(7) 


by  the  mean  value  theorem.  Over  this  region  then  (using  the  redefined  t)  the  equation  for 
the  foundation  velocity  becomes 


Z(T)  =  Zn 


S„t 
+  TT 


(8) 


♦This  is  an  unnecessary,  but  convenient  restriction.  It  makes  the  calculations  less 
cumbersome  and  machine  programming  easier. 
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e- 


where  sn  =  Zn  +  t  -  Z  ■  sn  Is  in  reality  a  first  forward  difference  and  Eq.  (8)  in  effect 
replaces  the  curved  foundation  velocity  function  in  this  region  by  a  straight  line.  If  greater 
precision  is  desired  an  assumption  can  be  made  about  the  rate  of  change  of  the  foundation 
acceleration,  such  as,  let  Z  be  constant  over  the  increment.  A  parabola  passed  through 
the  points  Zn-1,  Zn,  and  Zn+1,  or  one  passed  through  Zn,Zn  +  , ,  Zn  +  2 ,  will  satisfy  this 
requirement  over  the  increment.  Therefore 


or 


Sn  t 

=Zn+-ir  + 
s„  t 

Z(t)  =  zn  + 


(9) 

(10) 


where 


Sn-1  =  Sn  ■  ®n-l  “  ^n+l  "  2  +  ^n-1 

Sn  =  Sn  +  1  -  Sn  =  ^n  +  2  "  ^  ^ni  1  +  ' 

Note  that  S2  is  not  the  square  of  S  but  is  the  second  forward  difference.  In  both  equations 

when  t  =  h,  Z(h)  =  2n  +  1,  so  the  curve  is  not  in  this  respect  like  a  Taylor  or  Maclaurin 

series  expansion  which  has  been  truncated.  It  is  generally  more  convenient  to  use  Eq.  (9) 
in  the  numerical  problem,  because  after  the  solution  is  started  it  uses  a  rearward  point, 
and  no  difficulties  are  encountered  at  the  end  of  the  forcing  function. 


Fig,  Z  -  Portion  of  a  velocity  record 


Over  the  time  interval  from  Zn  to  Zn+1  the  equation  of  the  curve  has  now  been 
defined,  so  the  integration  and  evaluation  at  h  can  be  accomplished.  If  the  point  by  point 
numerical  integration  equations  are  then  scaled  by  w  for  convenience  in  number  handling, 
the  results  are 

.  S_  (1  -  cos  oj  h) 

*n  +  i“  =  Xnucoswh  +  Xnsin«h - — r - 


( 


1  t  cos  wh 


3  in  a  hN 
«2  h*  ) 


2  wh 


(Ha) 
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and 


Xn  +  1  =  -  Xnw  sir.  coh  +  Xn 


<uh 


Sn  s in  (o h 
o;h 


pj  ( 1  -  cos  ctih 

-1  l  *2h2 


sin  o>h 
2  o-h 


dib) 


Two  reasons  for  using  equally  spaced  time  increments  are  apparent  immediately:  (a)  it  is 
comparatively  easy  to  use  the  parabolic  approximation  of  the  forcing  function  with  equally 
spaced  points;  (b)  for  any  frequency  the  sines  and  cosines  need  only  be  calculated  once,  and 
hence  the  trigometeric  functions  merely  become  numerical  constants  multiplying  Xk,  X,  S, 
and  s2-  Of  course  if  the  parabolic  approximation  to  the  forcing  function  is  not  desired,  the 
32_j  terms  can  be  ignored,  or  higher  order  approximations  to  the  input  can  be  considered. 
The  natural  expansion  functions  for  the  response  of  an  oscillator  are  sines  and  cosines 
(or,  in  the  damped  case,  damped  sines  and  cosines)  and  the  results  are  of  this  form. 

There  still  remains  the  problem  of  starting  the  solution.  At  the  beginning  of  the  input 
Eqs.  (6)  are  still  valid.  It  is  only  necessary  to  account  for  the  initial  conditions,  if  they 
are  different  from  zero,  and  for  the  second  difference  term,  s| ,  which  does  not  exist 
for  the  first  step  of  the  solution.  Assume  the  general  initial  conditions  that  at  the  begin¬ 
ning  of  the  forcing  function  (in  reality  at  the  left-hand  limit  of  (t  -  0)  as  t  -  0,  that  is,  as 
t  -0  from  the  negative  direction)  that  x  =  d,  X  =  v,  and  there  is  a  step  change  in  foundation 
velocity  V.  Then  the  starting  equations  become 


and 


SD  (1  -  cos  oih) 

Xj  ai  -  do.  cos  u)h  +  (v  -  V)  sin  wh - - 


„2  /I  -  cos  t..h  sin  oih\ 

°  \  2“h  J*  h2  / 


(12a) 


X 


1 


=  -  di.  sin  «h  +  (v  -  V)  cos  «h  -  Su 


sin  c^h 

cJK 


2  (\  -  cos  u-h  sin  u,h\  , 

°  l  w2  h2  "  _2"h  /  ’ 


(12b) 


Since  during  the  integration  each  pair  of  points  x™  and  x  are  considered  to  be  new  “Initial 
conditions”  for  the  next  points,  these  equations  also  can  be  used  to  handle  finite  discon¬ 
tinuities  in  the  foundation  velocity,  at  any  time. 


The  step  by  step  numerical  integration  equations  exactly  satisfy  the  original  differ¬ 
ential  equation  over  any  increment  which  is  a  straight  line  (s2_j  made  zero)  or  a  parabolic 
arc.  Therefore,  consider  the  problem  of  the  response  of  a  linear  oscillator  to  a  velocity 
forcing  function  which  has  at  most  a  finite  number  of  finite  discontinuities  and  the  segments 
between  these  discontinuities  are  a  combination  of  straight  lines  and  second  order  parab¬ 
olas.  The  proper  use  of  these  numerical  integration  equations  will  exactly  solve  the 
problem  at  each  calculated  point  with  the  exception  of  number  round-off  error,  of  course. 
The  numerical  integration  equations  are  stable  and  exact  for  this  case  regardless  of 
increment  size.  For  segments  which  are  higher  order  curves,  the  differential  equation 
Is  exactly  satisfied  up  to  the  second  difference  in  the  forcing  function.  The  closer  the 
approximating  curve  lies  to  the  true  one,  the  more  precise  the  answer.  The  error  In 
this  technique  lies  in  the  closeness  of  fit  of  the  approximating  segmental  curves  to  the 
forcing  function;  no  error  is  Inherent  due  to  Increment  size  alone,  in  the  linear  case. 
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Damped  Response 

The  differential  equation  of  motion  to  be  solved  for  the  damped  case  is 

x  +  an#  x  +  w2x  =  -  z  (5) 

where  again  it  is  assumed  the  foundation  velocity  is  the  known  forcing  function.  This 
differential  equation  has  the  solution  (1) 

X  =  X(0)  e~acJt  (cos  pt  +  ^=a_2  sin  pt^ 

t 

+  2l  -i  f  Z(T)  e'a“  [•  "tl  sin  p(t  -  T)  dT  .  (13) 

o 

If  the  same  procedure  is  used  to  approximate  the  forcing  function  as  in  the  damped  case 
then  the  step  by  step  numerical  integration  equations  (scaled  by  o>  for  convenience)  become 


Xn  e~awh  sin  ph 


■=  X„uc'<lMh  /cos  ph  +  T==  sin  ph\  +  ~ - - 

l  /l-a2  )  Vi  -a2 

-  S  [*  -  e_a"h  (cos  ph  +yfrp  sin  ph)] 

%{!  -  S  [ft  *  &)  -  l) 


xn+1  =  -  ♦  *„  e‘owh  (cos  ph  -  sin  ph) 

Sn  e_a"hsinph  Sn-l  \  1  _-awh  l-os  ph  Ja  1^  sinph]\,(UM 

■  a;  “  ’  ~  FR  "  nar-  +  ]y  ’ 

The  procedure  for  obtaining  the  starting  equations  is  the  same  as  the  undamped  case. 

A  Word  of  Caution 

Often  it  is  assumed  that  the  damped  and  undamped  natural  frequencies  are  equal.  If  such  an 
assumption  is  made  the  second  difference  (s2^)  should  be  dropped  from  the  calculation  because 
the  error  in  the  exponential  and  trigonometric  coefficient  of  this  term  in  the  X  equation  can 
become  larger  than  the  true  value  of  the  coefficient.  The  abridged  equations  become 

xn  +  1  oi  =  Xnoie-aa,h  (cos  tub  +  a  sin  wh)  +  Xn  e~a“h  sin  aih 

Sn 

-  ^  {l  -  e~awh  (cosiih  +  a  sin  wh))  (15a) 

Xn  +  1  =  -  Xn  <*>  e_awh  sin  wh  +  Xn  e~  a“h  (cos  wh  -  a  s  in  &ih) 

Sn 

n  .-auh 
“  ajh  6 


sin  coh  . 


(15b) 
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SAMPLE  PROBLEMS 

Two  simple  problems  were  chosen  to  show  sample  computations.  Known  solutions 
exist,  so  that  comparison  with  theoretical  values  is  possible.  Consider  the  following 
velocity  forcing  function:  Let 


Z=0,  for 

t 

S  0 

Z  “  Zq  s in  \ 

t  , 

for 

Z  =  0,  for 

t 

£  4T7/\ 

and  assume  the  initial  conditions  are  such  that  the  mass  is  at  rest.  The  values  of  the 
chosen  parameters  are  a>  =  600  rad/sec,  \  =  1200  rad/sec,  h  =77/7200  sec,  fa  =  95.5  cps, 
f  =  191.0  cps,  Z„  =  90 in. /sec, and  a  =  0.1  =  10%  of  critical  (when  damping  present). 

Table  1  is  the  listing  of  the  velocity  forcing  function  for  24  increments  needed  for  the 
integration.  The  values  of  the  sines,  cosines,  and  exponentials  used  in  these  examples 
were  taken  to  six  decimal  places,  but  the  actual  numerical  integration  was  held  to  three 
decimal  places  in  the  answers. 


Table  1 

Velocity  Forcing  Function 


n 

Sn  _  _ 

c2 

an-l 

0 

0 

,+45.000  000 

-12.057  714* 

1 

+45.000  000 

+32.942  286 

-12.057  714 

2 

+77.942  286 

+12.057  714 

-20.884  572 

3 

+90.000  000 

-12.057  714 

-24.115  428 

4 

+77.942  286 

-32.942  286 

-20.884  572 

5 

+45.000  000 

-45.000  000 

-12.057  714 

6 

0 

-45.000  000 

0 

7 

-45.000  000 

-32.942  286 

+12.057  714 

8 

-77.942  286 

-12.057  714 

+20.884  572 

9 

-90.000  000 

+12.057  714 

+24.115  428 

10 

-77.942  286 

+32.942  286 

+20.884  572 

11 

-45.000  000 

+45.000  000 

+12.057  714 

12 

0 

+45.000  000 

0 

13 

+45.000  000 

+32.942  286 

-12.057  714 

14 

+77.942  286 

+12.057  714 

-20.884  572 

15 

+90.000  000 

-12.057  714 

-24.115  428 

16 

+77.942  286 

-32.942  286 

-20.884  572 

17 

+45.000  000 

-45.000  000 

-12.057  714 

18 

0 

-45.000  000 

0 

19 

-45.000  000 

-32.942  286 

+  12.057  714 

20 

-77.942  286 

-12.057  714 

+20.884  572 

21 

-90.000  000 

+12.057  714 

+24.115  428 

22 

-77.942  206 

+32.942  286 

+20.884  572 

23 

-45.000  000 

+45.000  000 

+  12.057  714 

24 

0 

♦This  is  really  Sq  , 
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Undamped  Response 

The  numerical  integration  equations  tor  the  undamped  response  become 

xn  +  ,  w  =  0,965926  Xn  a  +  0.258819  Xn  -  0.130154  Sn  +  0.021593  (16a) 

Xn  +  l  =  -0.258819  Xn  ♦  0.965926  Xn  -  0.988616  Sn  -  0.002843  S^_,  (16b) 

Table  2  was  prepared  and  the  integration  carried  out  in  it.  The  theoretical X« is  included 
to  show  the  precision  of  the  solution.  The  forcing  function  was  defined  for  every  30  degrees 
(on  its  own  sine  wave)  and  the  response  was  therefore  calculated  at  a  time  which  corre¬ 
sponds  to  an  angle  of  15  degrees. 


Tabic  2 

Undamped  Response  Computations 


71 

0.965926 

0.258819 

-0.130154  ' 

*0.021593 

r 

Theoretical 

-0.258819 

♦0.965926 

-0.988616 

-0.002843 

L.  *■.  . 

Sn 

S*-, 

X„. 

. 

Xn 

*n+l 

0 

0 

0 

-5.857 

-0.260 

-6.117 

-5.994 

0 

0 

-44.488 

>0.034 

-44.454 

i 

-5.909 

-11.506 

-4.288 

-0.260 

-21,963 

-21.962 

♦1.583 

-42.939 

-32.567 

♦  0.034 

-73.889 

2 

-21.215 

-19.124 

-1.569 

-0.451 

.42.359 

-42.42? 

*5.684 

-71.371 

-11.920 

♦0.059 

-77.548 

3 

-40.916 

-20.071 

♦1.569 

-0.521 

-59.939 

-60.000 

♦  10.963 

-74.906 

♦11.920 

♦0.069 

-51.954 

4 

-57.897 

-13.447 

♦  4.288 

-0.451 

.67.507 

.67.491 

*15.513 

-50.184 

♦32.567 

♦0.059 

-2.045 

5 

-65.207 

-0.529 

♦  5.857 

-0. 260 

•  60.139 

-60.000 

♦17.472 

-1.975 

♦44.488 

♦0.034 

♦G0.010 

6 

-58.090 

♦15.534 

♦  5.857 

0 

-36.699 

-36.433 

♦  15.565 

♦57.974 

+44.488 

0 

+118.02? 

7 

-35.449 

.30.548 

♦4.288 

*0.260 

-0.353 

0 

♦  9.498 

♦114.005 

♦  32.567 

-0.034 

♦156.036 

8 

-0.341 

♦  40.385 

♦1.569 

*0.451 

*42.064 

*42.427 

*0.091 

•150.719 

•11.920 

-0.059 

♦162.671 

9 

♦40.631 

♦42.102 

-1.569 

♦0.521 

♦81.685 

♦  81.962 

-10.887 

♦157.128 

-11.920 

.0.069 

♦134.252 

10 

♦78.902 

♦  34.747 

-4.288 

*0.461 

♦109.812 

.100.918 

-21.142 

•129.677 

-32.567 

-0.059 

+75.909 

11 

♦106.070 

♦19.647 

-5.857 

•  0.260 

♦120.120 

*120.000 

-28.421 

♦73.322 

-44.488 

-0.034 

+0.379 

12 

♦116.027 

♦  0.098 

-5.857 

0 

•110.268 

♦109.918 

-31.089 

♦0.366 

-44.488 

0 

-75.211 

13 

♦106.511 

-19.466 

-4.288 

-0.260 

♦82.497 

♦81.962 

-28.539 

-72.648 

-32.567 

♦  0.034 

-133.720 

14 

♦79.686 

-34.609 

-1.569 

-0.451 

♦43.057 

♦42.427 

-21.352 

-129.164 

-11.920 

.0.059 

.162.377 

15 

♦41,590 

-42.026 

♦  1.569 

-0.521 

♦  0.612 

0 

.11.144 

-156.844 

♦11.920 

♦0.069 

-155.999 

16 

♦0.591 

-40.378 

.4.288 

-0.451 

-35.948 

-38.433 

-0.158 

-150.683 

♦  32.567 

+0.059 

-118.215 

17 

-34.723 

-30.596 

♦  5.857 

-0.260 

.59.722 

-60.000 

♦  9.304 

-114.187 

♦44.488 

+0.034 

-60.361 

18 

-57.687 

-15.623 

♦  5.857 

0 

-67.453 

-67.491 

♦  15.457 

-58.304 

♦44.488 

0 

+1.641 

19 

-65.155 

♦0.425 

•4.288 

♦0.260 

.60.182 

-60.000 

♦  17.458 

♦1.585 

♦32.567 

-0.034 

+51.576 

20 

-58.131 

.13.349 

♦1.569 

♦  0.451 

.42.762 

-42.427 

♦  15.576 

*49.819 

♦  11.920 

-0.059 

+77.256 

21 

-41.305 

♦  19.935 

-1.569 

♦0.521 

-22.358 

-21.962 

*11.068 

♦  74.624 

-11.920 

-0.069 

+73.703 

22 

-21.596 

♦  19.076 

-4.288 

♦0.451 

-6.357 

.5.994 

♦  5.787 

.71.192 

-32.5G7 

-0.059 

+44.353 

23 

-0.140 

1  ♦11.479 

-5.857 

♦  0.260 

-0.258 

0 

♦  1.645 

♦42.842 

-44.488 

-0.034 

-0.035 

Damped  Response 

For  this  problem  the  numerical  equivalent  of  Eqs.  (14)  become 

xn  +  l  ca=  0.966512  Xn  u  +  0,252160  X„  -  0.127915  Sn  +  0.021036  (17a) 

Xn  +  1  =  -0.252160  Xn  to  +  0.916080  Xn  -  0.963180  Sn  -  0.007006  s£_j  (17b) 
The  results  and  comparisons  with  the  theoretical  solution  are  shown  in  Table  3. 


10 


NAVAL  RESEARCH  LABORATORY 


Table  3 

Damped  Response 


" 

0.966512 

0.252160 

*n 

-0.127815 
Sn _ 

0.021036 

JL i_ 

Theoretical 

-0.252160 

X„* 

0.916080 

xn 

-0.963100 

-0.007006 
.  S^l__ 

0 

0 

0 

-5.756 

-0.254 

■  6.010 

-5.800 

0 

0 

-43.343 

+0.084 

-43.259 

1 

-5.809 

-10.908 

-4.214 

.0.254 

-21.165 

-21.190 

♦  1.515 

-39.629 

-31.729 

♦0.084 

-69.759 

2 

-20,476 

-17.590 

-1.542 

-0.439 

-40.047 

-40.187 

♦5.342 

-63.905 

-11.614 

♦0.146 

-70.031 

3 

-38.706 

-17,659 

+1.542 

-0.507 

-55.330 

-55.301 

♦  10.098 

-64.154 

+11.614 

+0.169 

-42.273 

4 

-53.477 

-10.660 

+4.214 

.0.439 

-60.362 

-60.344 

♦  13.952 

-38.725 

+31.729 

+0.146 

+7.102 

5 

-58.341 

+1.791 

+5,756 

-0.254 

-51.048 

-50.909 

♦  15.221 

♦6.506 

+43.343 

♦0.084 

♦65.154 

6 

-49.339 

+  16.420 

+5.756 

0 

-27.154 

-26.896 

+12.872 

♦59.686 

♦43.343 

0 

+  115.901 

7 

-26.245 

-+29.226 

+4.213 

+0.254 

+7.448 

+7.779 

♦6.847 

.106,115 

+31.729 

-0.084 

♦  144.667 

8 

+7, 199 

+36.479 

+  1.542 

+0.439 

+45.659 

♦45.985 

-1.878 

♦132.527 

.11.614 

-0.146 

♦142.117 

9 

+44.130 

+35.836 

-1.542 

+0.507 

♦78.931 

♦79.183 

-11.513 

♦130.191 

-11.614 

-0.169 

+  106.895 

10 

+  78  288 

+26.954 

-4.214 

♦0.439 

♦99.467 

♦99.528 

-19.903 

♦97.924 

-31.729 

-0.146 

♦46.146 

11 

+96.136 

+11.636 

-5.756 

♦0.254 

♦  102.270 

♦102.119 

-25.082 

.42.213 

-43.343 

-0.084 

•  26.236 

12 

+98.845 

-6.616 

-5.756 

0 

♦  86.473 

♦86.117 

-25.788 

-24.034 

-43.343 

0 

-93.165 

13 

+83.577 

-23.492 

-4.214 

-0.254 

♦55.617 

♦55.113 

-21.805 

-85.347 

-31.729 

♦  0.084 

-138.797 

14 

+53.754 

-34.999 

-1.542 

-0.439 

♦  16.774 

+16.214 

-14.024 

-121.149 

-11.613 

♦  0.146 

-152.640 

15 

+16.212 

-38.490 

+  1.542 

-0.507 

-21.243 

-21.749 

-4.230 

-133.330 

.11.613 

♦0.169 

-132.278 

16 

-20.532 

-33.355 

+4.214 

-0.439 

-50.112 

-50.469 

♦5.357 

-121.177 

+31.720 

+0.146 

-83.945 

17 

-48.434 

-21.168 

+5.756 

-0.254 

-64.100 

-64.245 

♦  12.636 

-76.900 

+43.343 

+0.084 

-20.837 

18 

-81.953 

-5.254 

+  5.755 

0 

-61.451 

-61.373 

♦16.163 

-19.088 

+43.343 

0 

♦40.418 

19 

-59.393 

♦10.192 

♦4.214 

♦0.254 

-44.733 

-44.741 

♦  15.495 

♦37.026 

♦31.729 

-0.084 

♦84.166 

20 

-43.235 

♦  21,223 

♦1.542 

♦0.439 

-20.031 

-19.662 

+  11.280 

+77.103 

♦11.614 

*0.146 

+99.851 

21 

-19.360 

+25.178 

-1.542 

♦0.507 

♦4.783 

*5.159 

+5.051 

♦91.472 

-11.614 

-0.169 

+84.704 

22 

+4.623 

+21.368 

-4.214 

+0.439 

♦22.216 

♦25.509 

-1.208 

♦77.629 

-31.729 

-0.146 

♦44.548 

23 

+21.472 

+  11.233 

-5.756 

♦0.254 

♦27.203 

♦27.206 

-5.602 

♦40.810 

-43.343 

-0.084 

♦8.219 

Increment  Size 

The  size  of  increment  for  shock  spectrum  calculations  is  restricted  if  the  method 
is  used  without  modification.  Since  this  is  a  point  by  point  technique,  if  h  =  Atls  chosen 
too  large  for  the  high  frequencies,  it  may  turn  out  that  values  near  the  peak  are  missed. 
This  results  in  a  value  which  is  too  low.  In  the  Fourier  portion  of  the  calculations  (to 
be  discussed  later)  the  only  restriction  on  increment  size  is  that  the  arcs  approximating 
the  forcing  function  must  lie  close  to  the  true  curve,  because  the  method  must  compute 
the  end  point  anyway.  A  good  rule  of  thumb  is  to  note  if  a  set  of  trapezoids  at  a  chosen 
increment  size  would  be  a  good  approximation  to  the  curve. 

For  forcing  functions  which  are  a  combination  of  straight  lines,  are  parabolic  arcs, 
or  are  “sufficiently  smooth,*  the  numerical  integration  equations  themselves  can  be  used 
for  interpolation  purposes*  without  changing  the  approximation  to  the  forcing  function. 

The  analyst  should  be  careful  however  not  to  use  the  values  at  n  +  1/2  when  integrat¬ 
ing  between  n  +  1  and  n  +  2.  This  technique  should  be  particularly  useful  when  hunting 
for  maximum  values  and  has  no  justification  or  need  in  Fourier  spectrum  type  of 
calculations. 

The  NRL  digital  computer  NAREC  has  been  coded  to  use  this  technique  (without 
interpolation)  in  the  case  of  the  velocity  inputs  and  some  results  are  shown  later. 


♦Caution;  This  ia  not  a  normal  linear  or  quadratic  interpolation.  See  Appendix  A. 
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FOURIER  ANALYSIS 
Fourier  Integrals 

Let  a  function  f(t)  be  real  and  satisfy  the  Dirichlet  conditions  (10),  and  let 


■fee 

I 


|f(t)|dt 


exist.  Then 


1/2  [f(t  +  0)  +  f(t  -  0)]  = 


''I  f 


f(s)  cos  w(s  -  t)  d3dw. 


(18) 


(19) 


This  is  known  as  the  Fourier  Integral  Theorem  (10-12).  The  function  on  the  left  merely 
indicates  that  the  double  integral  converges  to  the  average  value  of  the  left-  and  right-hand 
limits  at  a  finite  discontinuity  of  f(t).  For  the  remainder  of  this  paper  this  will  be  written 
as  f(t) ,  but  this  point  should  be  remembered. 

If  f(t)  Is  assumed  to  be  either  an  even  function  (f(-t)  =  f(t))  or  an  odd  function 
(f(_t)  =  simplifications  are  possible.  They  result  in 


Ht) 


CD  tO 

ai 


f(s)  cos  his  cos  cot  dsdu 


for  f(t)  even  and 


f(t) 


ao  too 

■iff 


f (S)  sinus  sin  ut  dsdu 


for  f(S)  odd. 

It  is  convenient  to  have  the  previous  Fourier  sine  and  cosine  integrals  as  transform 
pairs.  The  cosine  transform  is, 


where 


and  the  sine  transform  is 


where 


Fc(«) 


f(s)  cos  us  ds 


f(t) 


<0 

■if 


Fc(a>)  cos  wt  daj 


?.(«)  = 

1  f  (s)  sin  a»s  ds 

0 

a> 

A 

f(t)  =  # 

F, (t<>)  sin  wt  dt 

(20a) 


(20b) 


(21a) 


(21b) 
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Examination  of  Eq.  (20b)  will  show  that  if  t  =  0 

00 

o 

This  fact  has  application  when  the  inverse  transform  of  a  function  is  to  be  obtained. 


Fourier  Series 

It  is  often  convenient  to  expand  an  arbitrary  function  in  the  range  from  t  =  0  to  t  =  T0 
in  a  trigonometeric  series  of  the  form 


.  2n  nt 

A„  cos  -ii —  + 


00 

L 


B_  sin 


2nmt 


(22) 


This  is  the  full  range  expansion  form  of  the  Fourier  series  and  the  coefficients  are  found 
by  means  of  the  formulas 


and 


f(t)  dt, 

f (t)  cos  dt 

1  o 


(23) 


(24) 


....  ■  2n  v\ 

f  ( t)  sin  -= — 
Jo 


dt. 


(25) 


Half  Range  Expansions 

However,  it  may  become  convenient  to  expand  in  a  series  of  only  sines  or  cosines. 
For  example 


or 


'«>■£ 


. ,  ,  A0  x- 1  —  n?rt 

f(t)=T+  2^  Anc»s  T7 


These  coefficients  are  given  by 


(26) 


(27) 


(28) 
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and 


To 

A0  =  ^j  f(t)dt 
0 


129) 


(30) 


Parseval’s  Theorems 

Of  the  many  theorems  associated  with  Fourier  analysis  some  of  the  most  useful  are 
known  as  Parseval’s  theorems.  These  are  presented  without  proof.  For  the  Fourier 
integrals, 


J  f2(t)dt  =fj  F*(^)  d»  =  fj  F*(.J>da.c  £  J  F(2>  da,. 


For  Fourier  series, 


To  To 

/•  t  A2  “ 


And  for  the  half  range  expansions, 
To 


/ Tf  -'it* 


FOURIER  ANALYSIS  AND  THE  LINEAR  OSCILLATOR 

The  integrals  of  the  Fourier  transform  pairs,  and  for  calculation  of  the  Fourier  series 
coefficients  are  all  of  the  form 


or 


[ 

L 

I 


g(r)  sin  nr  dr 


g(r)  cos  ar  dr . 


(31a) 


(31b) 


Now  sine  and  cosine  functions  are  obviously  the  natural  expansion  functions  when  dealing 
with  the  solution  of  an  undamped  linear  oscillator.  The  question  to  be  answered  is  then, 
“Is  it  possible  and  efficient  to  use  the  undamped  response  of  a  linear  oscillator  to  an 
applied  transient  to  calculate  Fourier  spectra,  Fourier  series  coefficients,  and  inverse 
Fourier  transforms?”  The  answer  is,  “Yes.” 
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Fourier  Transform 

To  examine  the  problem  in  detail  it  is  now  time  to  return  to  Eq.  (6a),  which  is  repeated 
for  convenience, 


Xn 

X  =  X0  cos  o,t  +  —  sin  cjt  -  —  I  Z(T)  sin  w(t  -  T)  d" 


"I* 


(8a) 


If  the  Duhamel  integral  is  integrated  by  parts  and  then  the  mass  is  assumed  at  rest 
initially,  there  results 


X  =  - 


j*  Z(T)  cos  w(t  -  T)  dT. 


Use  of  the  trigonometrical  substitution 


cos  (a  -b)  =  cos  a  cos  b  +  sin  a  sinb 


(32) 


yields 


X  =  -  cos  wt  J”  Z(T)  cos  wt  dT  -  s in  wt  j* 


Z(T)  sin  uT  dT. 


(33) 


Both  of  these  integrals  can  be  associated  with  the  Fourier  transform  pairs  for  time  equal 
to  t,  so  the  equation  may  be  written 


X  =  Vf^(cu)2  +  F,(<o)2  cos(«t-0) 


where 


0  =  tan-1 


£. 

Fr 


Fc 

i*\  +  *1 ' 


The  quantity  (F*  t  F* )'  '2  is  usually  called  the  Fourier  spectrum  magnitude  value  at  u, 
and  0  is  the  phase  angle  associated  with  this  vibration.  The  maximum  value  of  this 
equation,  evaluated  after  the  shock  motion  is  over,  yields  the  interesting  result  that 
the  “after"  or  “residual”  shock  spectrum  is  In  reality  the  Fourier  spectrum  of  the 
foundation  velocity. 


Equation  (33)  and  its  derivative  give 


X  =  -  cos  «Jt  J  Z(T)  cos  a<T  dT  -  s in  cut  j*  Z 


Z(T)  s  i  n  o>T  dT  . 


and 


Z(T)  cos  oiT  dT  -  cos  o)t  I  Z(T)  sin  wT  dT. 


(34a) 


(34b) 
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If  X,  X,  and  Z  are  known  for  a  particular  time  T0  (in  this  case,  say,  the  end  of  the  transient), 
then  simultaneous  solution  of  these  equations  gives  the  integral  values.  They  are 

•  * 

j  Z(T)  cos  MT  dT  =  -  X(T0)  cos  mT0  +  X(T0)  +  Z,(T°^  sin  tiT0  (35a) 

Jo 

and 

•  X(Tn)  +  Z(Tn) 

Z(T)sinwTdT  = - ^ - -  cos  a>T0  -  X(T0)  sin  uT0  .  (35b) 

However,  these  are  precisely  Fc(<*>)  and  F,(<^)  as  defined  by  Eqs.  (20)  and  (21),  for  a  func¬ 
tion  which  ends  in  finite  time  and  are  of  the  form  of  Eqs.  (31). 


Fourier  Series  Coefficients 

To  find  Fourier  series  coefficients  for  the  full  range  expansion  let  *>„  =  2nw/T0,  and 
evaluate  the  response  X  and  X  at  T0  for  each  n  desired. 


Then 


and 


.  _  2  I  u,_. _ 2n-m  2X 

An  -T“  Z(T)COS^= - dT^-Tp- 

‘0  J  *0  l0 


To 

Bn=^j 


Z(T)  sin^Llil  dT=  -■£- 

T0  T0  “n 


For  the  half  range  expansions  let  on  =  nw/T0.  The  coefficients  become 


and 


An  =-  /  (-1)"  X 


5  2_  ,  ,n  (X  +  Z) 

“n  "  T-  '  *)  — 77, - 


where  the  x  and  x  have  again  been  evaluated  for  each  n  at  t  =  T0 .  This  is  an  automatic 
consequence  of  Eqs.  (35). 


Inverse  Fourier  Transform 

With  the  exception  of  a  constant  multiplier,  the  transform  and  its  inverse  as  defined 
here  are  symmetric.  Therefore,  if  one  is  willing  to  stop  the  &>  integration  short  of 
infinity,  an  approximation  to  the  inverse  transform  can  be  found  by  the  same  methods. 
However  some  unusual  scaling  is  involved,  especially  if  a  machine  program  has  been 
used  which  scales  the  x  and  x/w  by  w.  When  the  direct  transform  of  an  arbitrary  function 
is  found,  Fourier  sine  and  cosine  spectra  are  scaled  by  o>  with  the  proposed  numerical 
integration  equations.  (This  was  so  that  more  convenient  numbers  could  be  handled  but 
is  not  a  restriction  as  they  could  be  de-scaled.)  Consider  the  Fourier  cosine  transform 
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only,  because  the  same  approach  applies  to  the  sine  transform.  Suppose  a  machine  program 
calculated  o>fc  as  a  function  of  frequency  f .  However,  for  the  input  let  g(T)  in  reality  be 
g(<n)  and  attempt  to  obtain  the  inverse  transform.  The  machine  at  each  frequency  f  calcu¬ 
lates 


L 

2n  f  Hc  =  2v  f  g(w)  cos  2^f  <o  dw 

•b 

whereas 

L 

g(aj)  cos  cot  dco 


is  desired.  If  the  integral  was  multiplied  by  1  /■« 2 f  for  each  f ,  then  the  coefficient  of  the 
integral  would  be  correct.  Time  t  in  the  second  equation  must  correspond  to  2nf  in  the 
first,  so  t  =  2-ni.  An  example  showing  the  application  of  the  numerical  integration  equa¬ 
tions  to  the  inverse  transform  problem  is  shown  later. 

There  are  of  course  many  examples  of  the  use  of  Fourier  transforms  in  the  literature. 
Two  of  the  interesting  ones  are  the  calculation  of  frequency  response  knowing  the  time 
response  to  impulse,  and  the  inverse  of  this,  namely,  the  calculation  of  the  time  response 
for  impulse  knowing  the  frequency  (steady -state)  characteristics.  To  demonstrate  this 
consider  the  Duhamel  integral 


g(w)  =  -  f  (T)  cos  a>  (  t  -  T)  dT. 

*b 

By  successive  changes  of  variable,  this  may  be  written  as 

t 

g{w)  =  -  f(t  -  T)  cos  <oT  dT. 

*1 

If  f  (T)  was  the  response  to  impulse  the  cosine  term  could  be  considered  to  be  the  driving 
function  and  the  response  to  steady -state  vibration  can  be  found  for  large  t. 


APPLICATION  TO  NONLINEAR  PROBLEMS 

The  numerical  integration  equations  as  derived  here  can  be  applied  to  nonlinear  single 
degree  of  freedom  systems  and  to  two  degree  of  freedom  linear  and  nonlinear  systems. 
This  has  been  done  and  will  be  the  subject  of  a  later  NRL  report.  To  merely  indicate 
the  method  of  solution  consider  a  nonlinear  oscillator  responding  to  an  applied  force.  The 
equation  of  motion  might  be  written  as 

X  +  2a wX  +  u>2X  =  g(t)  -D(X)  -  K(X).  (2a) 

Now  if  over  a  short  time  increment  the  quantity  on  the  right-hand  side  was  considered  to 
be  -An  which  is  defined  as  the  negative  of  the  average  value  over  this  Increment  then 
this  equation  may  be  written  as 
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This  over  this  short  time  interval  Is  a  second -order  linear  differential  equation  with 
constant  coefficients  and  the  solution  is  well  known.  Since  the  average  value  of  An  is 
required,  iteration  can  be  used  with  the  proper  numerical  integration  equations  as  derived 
here,  and  a  step  by  step  evaluation  of  the  response  is  possible. 

The  undamped  response  equations  could  also  be  used  by  putting  the  x  term  on  the 
right-hand  side.  This  then  becomes  the  numerical  equivalent  of  the  phase  plane  graphical 
method . 


EXAMPLES  OF  MACHINE  COMPUTATIONS 
Fourier  Series  Coefficients 

The  20th  to  the  30th  Fourier  cosine  series  (half range  expansion)  coefficients  of  the 
function  y  =  («/4)  I  sin  2^t  |  are  used  to  illustrate  this  application.  The  arch  of  the  sine 
curve  was  divided  into  4,  6,  8,  10,  and  12  increments  (increment  sizes  1/8,  1/12,  1/16, 

1/20,  and  1/24  the  size  of  a  sine  period)  and  for  each  of  these  increment  sizes  the  response 
of  an  undamped  linear  oscillator  was  used  to  calculate  the  coefficients  at  2-cps  increments 
from  40  to  60  cps.  The  input  was  given  to  six  significant  figures.  Table  4  shows  how  many 
periods  of  the  responding  system  correspond  to  the  increment  size  used.  For  example, 
at  60  cps  when  10  increments  were  used  (n  =  10)  the  ratio  of  the  increment  size  (1/20)  to 
the  period  of  the  responding  oscillator  (1/60)  is  3.0. 


Table  4 

Number  of  Periods  of  the  Responding  System 
Corresponding  to  the  Increment  Size  Used 


Oscillator 

Frequency 

(cps) 

Ratio  of  Increment  Size  to  Oscillator  Period 
When  the  Number  of  Increments  Used  is  n 

n  =  4 

n  =  6 

n  =  8 

n  =10 

n  =12 

40 

5.000 

3.333 

2.500 

2.000 

1.167 

42 

5.250 

3.500 

2.625 

2.100 

1.750 

44 

5.500 

3.667 

2.750 

2.200 

1.833 

48 

5.750 

3.833 

2.875 

1.917 

48 

6.000 

4,000 

3.000 

2.400 

2.000 

50 

6.250 

4.167 

3.125 

2.500 

2.083 

52 

6.500 

4,333 

3,250 

2.600 

2.167 

54 

6.750 

4.500 

3.375 

2.700 

2.225 

56 

7.000 

4.667 

3.500 

2.800 

2.333 

58 

7,250 

4.833 

3.625 

2.900 

2.417 

7.500 

5.000 

3.750 

3.000 

2.500 

Table  5  contains  the  results  of  the  computation.  The  values  of  each  of  the  coefficients 
have  been  multiplied  by  minus  106.  The  numbers  have  been  rounded  to  the  sixth  decimal 
place  to  correspond  to  the  six  decimal  place  input.  To  show  the  precision  of  the  method 
the  errors  in  percent  of  the  calculated  value  when  compared  to  the  exact  theoretical  value 
are  given.  For  example,  at  60  cps,  with  n  =  10,  a30  =  0.000287  (calculated),  a30  =  0.000278 
(exact  theory),  and  Error  =  3,2%. 
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This  was  using  an  increment  size  which  is  three  times  the  responding  oscillator’s 
period,  Therefore  thirty  or  more  coefficients  can  be  calculated  using  only  a  few  input 
points. 

For  those  forcing  functions  which  are  exactly  satisfied  by  the  approximations,  the 
only  error  present  is  the  normal  round-off  error. 


Tabic  5 

Comparison  of  Calculated  and  Theoretical  Values  (Multiplied  by  10M  for  Fourier  Coefficients 


Frequency 

Theoretical 

Computed 

Percent 

Computed 

Percent 

Computed 

Percent 

Computed 

Percent 

Computed 

Percent 

U*ps) 

Value 

=  4 

Error 

n  =  6 

Error 

n  =  8 

Error 

n  =  10 

Error 

n  =  12 

Error 

40 

625 

728 

16.5 

660 

5.6 

663 

6.1 

645 

3.2 

634 

1.4 

4  2 

567 

693 

22.2 

580 

2.3 

581 

2.5 

608 

7.2 

580 

2.3 

44 

517 

544 

5.2 

547 

5.8 

543 

5.0 

539 

4.3 

532 

2.9 

46 

473 

574 

21.4 

534 

12.9 

509 

7.6 

484 

2.3 

487 

3.0 

48 

434 

505 

16.4 

470 

8.3 

455 

4.8 

43S 

i.i 

444 

2.3 

50 

400 

489 

22.3 

458 

14.5 

439 

9.8 

403 

0.8 

422 

5.5 

52 

370 

390 

5,4 

391 

5.7 

389 

5.1 

375 

1.4 

383 

3.5 

54 

f  43 

417 

21,6 

351 

2.3 

35) 

2.3 

352 

2.6 

351 

2.3 

56 

319 

371 

16.3 

337 

5.6 

323 

1.3 

332 

4.1 

323 

1.3 

58 

297 

363 

22.2 

336 

13.1 

305 

2.7 

312 

5.1 

300 

1.0 

60 

278 

293 

301 

8.3 

292 

5.3 

287 

3.2 

279 

0.4 

Inverse  Fourier  Transforms 

To  demonstrate  that  this  method  can  be  used  to  calculaie  the  inverse  Fourier  trans 
form  without  the  use  of  special  tables  the  following  two  cases  were  worked  out  on  the 
digital  computer  NAREC. 


Case  I  -  The  Fourier  transform  of  a  step  function  of  heigth  one-half  unit,  lasting  two 
seconds  of  time,  and  returning  to  zero  with  a  step  change  of  one-half  unit  was  used  as  an 
input  to  the  shock  record  program.  (This  program  prints  out  the  scaled  Fourier  transform, 
the  Fourier  sine  and  cosine  transform,  as  well  as  +  Xw(max)  and  -Xo>(max).)  Of  course  the 
integration  could  not  be  carried  out  to  infinite  frequency  so  a  cutoff  frequency  of  approxi¬ 
mately  83  cps  was  decided  upon  and  about  100  output  points  were  used  to  define  the 
function,  up  to  time  equal  to  v  seconds.  The  error  in  the  resulting  transformation  was  so 
small  that  a  graphical  plot  of  the  original  function  superposed  upon  the  results  of  the 
tnversc  transformation  will  not  show  the  precision  of  the  method.  Table  6  shows  the 
average  algebraic  error,  the  average  absolute  error,  and  the  root  mean  square  error. 

Case  II  -  As  possibly  a  more  severe  test  case,  it  was  decided  to  try  a  damped  sine 
wave.  The  procedure  was  to  let  the  NAREC  calculate  the  original  function  for  about  four 
complete  cycles,  defining  it  at  100  points.  This  digital  record  was  then  used  as  an  input 
to  the  shock  record  program  to  calculate  the  Fourier  sine  and  cosine  transforms  at 
intervals  ol  2  cps  from  0  cps  to  300  cps.  The  sine  and  cosine  transforms  were  then  used 
to  calculate  the  original  time  function  again  by  using  them  as  inputs  to  the  shock  record 
program,  The  original  time  function  was 

y  =  e"a<lJt  sin  pt 

where  a  =  0.2,  u  =  300,  (f  ~  47  cps),  p  =  a>/i  -  =  293.938+,  and  the  maximum  value  of 

y  turned  out  to  be  about  0.738.  Figure  3  shows  the  results  of  this  computation.  Agreement 
would  have  been  better  if  a  higher  cutoff  frequency  had  been  chosen,  say  500  cps  or  750 
cps,  and  if  more  closely  spaced  points  had  been  used  to  define  the  original  Input  function. 
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Table  6 

Numerical  Results  for  Case  I 


Error 

Computation* 

Inverse  Fourier  Cosine 
Transformation 

Inverse  Fourier  Sine 

Transformation 

... 

Average  _  21  d ; 
Error  *N 

0.23  x  10'3 

0.14  x  10~3 

Average  £|d.| 

Absolute  =  j _ 

Error  N 

0.95  x  10‘3 

1,33  x  10'3 

Hoot  Mean  /=~ 
Square  =  d  i 

Error  N 

0.1B  x  10-3 

0.28  x  10'3 

*dj  is  the  deviation  at  point  i  and  N  is  the  total  number  of  data 
points  (101  for  the  cosine  transform,  and  94  for  the  sine 
transform), 


Fig.  3  -  Comparison  of  results  for  Case  II,  a 
damped  sine  wave 
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Damped  Shock  Spectra 

To  show  the  effect  of  damping  upon  shock  spectra  an  actual  foundation  velocity  record 
(Fig,  4)  was  used.  It  was  divided  into  760  increments  and  shock  spectra  for  damping  ratios 
of  a  =  0,  a  =  0,01,  and  a  =  0.10  were  calculated  and  are  shown  in  Fig.  5.  (Curves  for  a  = 
0.001  and  a  =  0.0001  were  also  obtained  but  are  not  shown  for  reasons  of  clarity.)  Figure  6 
is  a  plot  of  the  effect  of  damping  upon  the  peak  value  and  the  effect  upon  the  shock  spectrum 
value  useful  in  stress  analysis  (4)  (that  value  corresponding  to  that  fixed  base  natural 
frequency  of  the  structure). 


Fig.  5  -  Undamped  and  damped  shock  spectra 
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Fig.  6  -  Effect  of  damping  upon  the  shock  spectrum 
values  of  interest 


SUMMARY  AND  CONCLUSIONS 

An  easily  understood,  yet  precise,  numerical  integration  method  has  been  derived  which 
allows  not  only  shock  calculations  to  be  made  but  also  Fourier  analysis.  This  technique, 
like  all  others,  has  two  kinds  of  error  present:  Inherent  and  round  off.  For  initial  value 
problems,  and  problems  with  forcing  functions  which  have  a  finite  number  of  finite  dis¬ 
continuities  and  can  be  exactly  described  by  a  set  of  straight  lines  or  parabolic  arcs,  the 
method  (as  handled  by  a  desk  operator)  has  no  Inherent  error.  For  continuous  curves 
the  differential  equation  is  exactly  satisfied  up  to  the  second  difference  of  the  forcing 
function.  The  Inherent  error  lies  in  the  closeness  of  fit  rather  than  in  any  built-in  ones. 
This  method  has  its  own  advantages  and  disadvantages  which  might  be  partially  listed  as 
follows. 

1.  High  precision  when  using  second  differences  since  error  depends  upon  “closeness” 
of  fit. 

2.  Response  including  a  linear  damping  assumption  is  not  much  more  difficult  to 
compute  than  response  with  no  damping  because  the  exponentials,  sines,  and  cosines 
must  be  calculated  only  once  per  frequency,  if  the  interval  is  kept  constant. 

3.  There  are  no  pseudoequations  assumed  for  the  original  differential  equation. 

4.  There  are  no  iterations  required  for  linear  systems. 

5.  The  numerical  solution  for  a  theoretically  solvable  problem  is  more  rapid  than 
the  exact  solution  because  the  original  values  of  the  trigonometric  and  exponential  func¬ 
tions  are  used  throughout  rather  than  being  calculated  anew  for  each  point. 

8.  As  has  been  indicated  and  will  be  shown  in  a  future  report,  the  method  can  be 
applied  to  nonlinear  one  and  two  degree  of  freedom  systems. 
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7.  As  w;th  am>  step  by  step  method  it  tends  to  become  tedious  for  a  lar~ge  number  of 
ii.crements  when  ..sing  a  desk  calculator  but  is  conveniently  coded  for  an  el  ectronic 
computer. 

8.  Each  new  value  depends  on  a  pair  of  previously  calculated  ones,  so  mistakes  are 
carried  -lone. 
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APPENDIX  A 


OTHER  FORMS  OF  NUMERICAL  INTEGRATION  EQUATIONS: 

INTERPOLATION 

The  numerical  integration  equations  are  presented  in  this  appendix  in  their  time 
form.  That  is  to  say  they  have  not  been  evaluated  for  t  =  h,  so  that  they  may  be  used 
to  interpolate  for  t  <  h.  To  obtain  the  normal  step  by  step  equations  let  t  =  h  and 
simplify.  The  input  (forcing  function)  is  defined  in  the  same  manner  as  in  the  derivation. 


FOUNDATION  VELOCITY  FORCING  FUNCTION 
The  undamped  equations  are 


Xn  +  1  "=  XnMCOS  cot  +  xn  sin  cot  -  Sn  (1  ~  ^  ~ 
c2 

n-1  /_t  ^  cos  cot  sin  e«>t  \ 

coh  \h  ”  2  +  2  roK  / 


and 


Xn  tl  =»XI1usinwt  +  Xn  cos  cot  - 


Sn  s in  cot 

coH 


Sn  - 1  (\ 

coh 


f  X.  cos  cot  sin  cot\ 
\coh  coh  2  1' 


The  damped  equations  are 


xn  +  i  "  =  Xn <o  e~a<0  *  (cos  pt  +-t=  sin  pt]  + 


Sj. 

coh 


Xn  e~aa>t  sin  pt 


and 


'n  +  l 


(  Vl  -  a2  J  yi  -  a2 

Tl  -  e~a“  1  /cos  pt  + - - —  sin  pt\ 

L  \  JJ 

^n-l  ft  1  2a  fa  ( 1  -  2a 2 )~]  e~aM*  sin  pt 
“coh  |h  “  2  ”  coh  +  [2  '  31  J  ^ - - 

+  (i  +  Ik)  e'a‘iJ,  cospt} 

/cos  pt  sin  pt\ 

V  i  i  -«*  ) 


-_ae-a<ot  sinpt 


ph 


coh 


cos  pt 


i 
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APPLIED  FORCES 

From  the  nth  to  the  (n  +  i)th increment  F(T)  is  defined  as 


F(T)  =  Fn 


S„T  s2_j 
+  — +  “ 


The  undamped  equations  are 

Xn  Fn 

Xn  +  i  =  Xn  cos  tut  +  —  sin  tut  +  -jj-  (1  -  cos  tut) 


and 


£n  /  t_  sin  tut\ 

+  K  \h  ‘  tuh  ) 

^n-1  ft2  t_  2(1-  cos  tut)  sin  cut 
+  2K  [h2  ‘  h  ‘  ffi2h2  +  tuh 


=  - Xn  co  s i n  cot  +  Xn  cos  cot  +  -j?  sin 

/l_**_cos_wt_\ 

*  K  \  <uh  ) 

Sn-l  f  2t  1  cos  tut  2  sin  tut  \ 

+  2K  ^2  ‘tuh +  ‘  tuh  ‘  w2h2  )  " 


The  damped  equations  are 


X  r  X  o_aw' 
An  +  1  -  An  e 


/cos  pt  +  a  —  sin  pt\  sin  pt 

\  ) 

-£r  jl  -  r.  ~a  "  ’  ^co  s  p  t  +  ^  a  —  sin  ptj 


Tf [f- Js  0 

Sn-lft2  t  4  a  t  [2(1  -4a2)  2a],, 

2K  |h2  cuh2  ’  L  to2  h  2 


cos  pt) 


fl  -  2a2  2a(3-4a2) 

e-QW*  sin  pt] 

+  tu2h2  J 

j 
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and 


y  V  p-aoit 

An+1  *n  e 


sin  pt  X„ 

n  0-ao;t 


^  1  -  a3 

Frt  -.-a  u)  t 


(cos  pt  -*■  a - sin  pt\ 

vrr^  ) 


"n  e'aMt  sin  pt  fl  e-a“*  /  a  .  \1 

«  •Kp'^r  ^r5""  jj 

pa  -  2d3)  _a_l  e~aMt  sin  pt~) 

‘  L  oj3  h3  "whJ  y,  _  a2~  ;• 


FOUNDATION  ACCELERATION 

By  means  of  the  similarity  between  Eqs.  (5a, b)  and  (5c, d)  it  is  only  necessary  to 
change  the  following  in  the  numerical  integration  equations  for  applied  forces  to  obtain 
those  for  foundation  acceleration,  Let 


F  =  -  M  Z„ 


K 


K 


and 


«$2 

sn-l 


APPENDIX  B 


ANOTHER  SET  OF  EQUATIONS  FOR  FOURIER  ANALYSIS 


The  solution  of  the  differential  equation  with  the  mass  initially  at  rest,  and  the  founda¬ 
tion  acceleration  as  the  driving  function  can  also  be  used  to  calculate  the  Fourier  Integrals. 
Consider 


and 


t 

X  =  -  i  Z(T)  sin  cu(t  -  T)  dT 
0 
t 

X  =-  J  Z(T)  cos  o>(t  -T)  dT. 


These  may  be  expanded  to 


and 


-  s  i  n  cut  f  Z 

;t  j  Z 


Z(T)  cos  mT  HT  +  cos  fi» t 


t 

J  Z(T)  sin  i 


>T  dT 


X  =  -  cos  uiK.  |  Z(T)  cos  wT  dT  -  s  in  cut  I  Z(T)  sin  o/T  dT. 

K 


t 

in  cut  J* 


Simultaneous  solution  yields 


f  z 


J 


Z(T)  cos  cuT  HT  =  -  Xcu  s  i  n  cut  -  X  cos  cut 


and 


Is 


Z(T)  sin  cuT  dT  =  Xcu  cos  cut  -  X  sin  cut . 
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